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SUMMARY 


A substructuring technique, originally developed for the effi- 
cient reanalysis of structures, is incorporated into the methodology 
associated with the plastic analysis of structures. An existing 
finite-element computer program that accounts for elastic-plastic 
material behavior under cyclic loading was modified to account for 
changing kinematic constraint conditions - crack growth and inter- 
mittent contact of crack surfaces in two dimensional regions. Ap- 
plication of the analysis is presented for a problem of a center- 
crack panel to demonstrate the efficiency and accuracy of the 
technique. 
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1. INTRODUCTION 


Fatigue crack propagation, until recently, was generally as- 
sumed to be directly related to the linear elastic stress-intensity 
factor range during cyclic loading. Implicit in this concept were 
the assumptions that only the tensile portion of the load cycle was 
effective in growing the crack and that crack closure occurs only 
at zero load. Elber (Refs. 1 and 2) has shown experimentally that 
fatigue cracks close at positive stresses during constant amplitude 
stress cycling. He has also indicated that fatigue-crack closure 
maybe a significant factor in causing the stress-interaction ef- 
fects (crack growth delay or acceleration) on crack growth under 
general cyclic loading. The crack closure phenomenon, associated 
with an extending crack, is believed to be caused by residual 
plastic deformations remaining in the wake of an advancing crack 
tip. A reasonable analytic model for crack closure and for the 
extending crack problem must possess the capability of accounting 
for changing boundary conditions (crack growth and intermittent 
contact of the crack surfaces) during a specified load history. 
These changing boundary conditions must be incorporated into the 
equations that govern the nonlinear load-deformation behavior. 

The present report is concerned with the modification of an 
existing nonlinear finite element program (Ref. 3) to account for 
crack extension and crack closure and the application of this pro- 
gram to the study of fatigue-crack closure. The modification of 
the existing finite element program consists of incorporating a 
method for the efficient reanalysis of structures (Refs. 4 and 5) 
that undergo changes in material properties or restraint condi- 
tions. The procedure has the advantage of hot requiring the re- 
formulation of the stiffness influence coefficients of the origi- 
nal, unmodified structure. 

The technique to treat plasticity in the FAST (Fracture 

Analys is~d3E~~STruc tures3~progr am-is—bas ed-on-the-init.ial._s t rain 

concept where an effective plastic load vector is introduced, in 
addition to the applied mechanical load, to account for the devel- 
opment of plastic deformation. Thus, the procedure used for the 
reanalysis of structures with variable restraint conditions is 
extended to include the effects of plasticity and has been incor- 
porated into the final program. 



Previous Studies 


The problem of an extending crack has been previously treated 
within the framework of a nonlinear finite element analysis. The 
procedure as described in Refs. 6 and 7 involves establishing an 
elastic-plastic state including the nodal forces (reactive forces) 
that hold together an element node directly ahead of a crack tip. 
An equal and opposite nodal force is then applied in increments to 
the crack tip node until this restraining force is completely re- 
moved. The node is then assumed to be free and displaces an 
amount representing the crack opening displacement. In Ref. 6, 
the crack tip is assumed to advance to a new position, correspond- 
ing to the adjacent node of the previous location of the crack 
tip. After a "step-of-growth, 11 as described in Ref. 7, the finite 
element mesh is shifted to a new position so that the new crack 
tip is in the same position as the previous one. Details asso- 
ciated with this shifting procedure, which simulates an infinite 
process, are not given. However, the following statement that 
appears in Ref. 6 should be noted: "In practice, the cumulative 

numerical errors involved in these incremental loadings become 
intolerable after a few small increments of crack extension." For 
reasons discussed in Section 2, the procedure used in Refs. 6 and 
7 appears to be incomplete and should, consequently, lead to 
numerical difficulties. 

Another approach (Ref. 8) , involves the use of a layer of 
"variable-material" elements beneath an axis of symmetry along 
which a crack can extend. These special elements take on soft, 
"jell-type" stiffness properties along the open (and extended) 
crack length and become relatively stiff, to simulate a rigid con- 
straint, during crack closure. This model, though physically sat- 
isfying, requires, at best a partial reformulation of the assem- 
bled stiffness matrix, and at worst (if the program is so con- 
structed) a complete reformulation and reanalysis. 

Another approach involves an iterative solution to the gov- 
erning equations that are determined for an a priori set of 
kinematic boundary conditions. The appeal of this approach is 
that the stiffness influence coefficient matrix remains unchanged. 
Thus, a solution technique such as the Cholesky triangularization 
of a positive definite, symmetric matrix need be performed only 
once, leading to a potentially rapid iterative process. However, 
application of this technique to simple yet representative crack 
extension problems reveals this procedure to be nonconvergent . It 
is our opinion that the process may be made convergent by the in- 
troduction of a relaxation factor in the iterative process. Since 
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the form of this factor is arbitrary and will, in general, vary 
with each problem, and within a problem for various combinations 
of separation and contact, the iterative approach was abandoned. 
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2. METHOD OF ANALYSIS 


The approach selected for the plastic analysis of structures, 
as available in the FAST program, incorporates the initial strain 
concept within the framework of the finite element method. In 
this technique the load-deflection relations are written to in- 
clude the effects of initial (or thermal) strains. These initial 
strains are then interpreted as plastic strains and the problem is 
solved by using the load-deflect ion- initial strain relations with 
subsidiary constitutive relations for an elastic-plastic engineer- 
ing material. An extensive presentation of the governing matrix 
equations and of the nonlinear constitutive relations is presented 
in Refs. 3, 9, and 10 and is only briefly reviewed here in the 
context of the modifications required to treat intermittent con- 
tact and separation. 


Governing Equations 

The matrix equation governing the response of a structure to 
some arbitrary history of loading can be written (as in Ref. 10); 


where 


[K]juj = jl 


+ |q 


+ |r! 


( 1 ) 


[K] = the conventional elastic stiffness influence 
coefficients 

{u} = generalized nodal displacements 
{P} = generalized nodal forces 


(Q) = "effective” plastic load that accounts for 
the presence and development of plastic 
strains 

{R} = an equilibrium imbalance that may exist as 
a result of the nature of the solution pro- 
cedure (Ref. 10) 

If we decompose the total strains, (e}, into elastic (e } 
and plastic components {eP} as 


4 


( 2 ) 



then for small deformations (no geometric nonlinearities), the 
stiffness matrix [K ] in Eq. (1) is the assemblage of the ele- 
ment stiffness matrices (see Ref. 10), and the effective plastic 
load vector {Q} is 


M 


i tk i ] ^ eP 

i=l 


r i 


(3) 


where _N is the number of finite elements in the plastic range, 
and [k] is the initial strain matrix of the individual ele- 
ments . Note that in the FAST program the plastic strains are 
assumed to be constant quantities (at centroids of elements) . 

If a distribution of plastic strain within an element is de- 
sired, its^ assumed functional form must be considered in deter- 
mining [k] . 

Solution Procedure for Unmodified Structure 

The algorithm for the incremental procedure used to solve 
for stresses, strains, and displacements for a typical load 
history follows: 

1) Determine the generalized displacements, {u} 
from the solution of Eq. (1) 

2) Use the solution from step 1 and appropriate 
strain-displacements relations to determine 
the total strains, {e}, and the increment 
of total strain [Ae] for the increment of 

load, {Ap} _ 

3) Compute the increments of stress {Ac} and 
the increment of plastic strain {Ae?} from 
the total strain increment {Ae} 

4) Determine an updated plastic load vector {Q} 
from Eq. (3) and mechanical load {p} k = 

{p} k " 1 + {Ap} 

5) Repeat steps 1 through 4 
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From the above it is apparent that at load step k the 
plastic strains and hence the effective plastic load vector 
cannot be determined without an a priori knowledge of the 
displacement. In the method chosen for use, the effective 
plastic load is taken to be equal to that computed in the pre- 
ceding load increment and is thus taken as a known quantity in 
the equation. 


The final form of Eq. (1) is 


[K]ju 


f I 1 

+ w 


where k and k-1 are the current and preceding load steps, 
respectively. The magnitude of the vector^ of equilibrium im- 
balance loads, (R) is a result of the incremental lineariza- 
tion of the nonlinear problem and the error associated with 
using an estimated value of [Q] . The effect of the lineariza- 
tion error is reduced by using sufficiently small load steps, 
whereas the latter error is accounted for by setting 


r} - X [k. ]| 


k. ]iAe p 


where {Ae P } is the increment of plastic strain determined in 
the preceding step. 

Solution Procedure for Variable Restraint Condition s 

^ It is desired to determine a set of modified displacements 
{u } = {u + 6u) due to a change in stiffness from [K] to 
[K + 6K], This modification may result from the separation 
(crack extension) or coalescing of two nodes to simulate con- 
tact. Equation (4) now becomes 


[K + 6K] ju* ■ = |f 


\ 


The modified displacement vector (u ) can, of course, be 
found by a complete reanalysis (including reassembly) using the 
modified stiffness influence coefficients. For localized 
changes (under consideration here) this would be most uneconomi- 
cal. 

The procedure chosen for implementation originally appeared 
in Refs. 4 and 5 and resembles the sub structuring approach pre- 
sented in Ref. 11. This technique, modified to treat the intro- 
duction of additional degrees of freedom, is incorporated into 
the FAST program, and is briefly described as follows: 

Consider Eq. (5) to be partitioned into 



where using the notation of Ref. 5, the subscripts i refer to 
the unmodified degrees of freedom and e to those that are to 
be modified; i.e., originally restrained against motion and now 
allowed to displace. Equation (7) may be expanded into the two 
equations : 

M u *l + M U e} -- M (8a) 


[K . ]ju*l 
er ) rj 


+ ' [K ]ju 
ee [ 


*' 

» 

e 



(8b) 


If we initially consider {u e j = 0, then we can, from 
Eq. (8a), solve for a quantity {tTjJ , as 

fi} = [K ii ] " I { F ij < 9 > 

\ . _ ^ 

By means of the Cholesky triangularization (see Appendix C) of 
a positive definite, symmetric matrix the stiffness matrix 
can be written as the product 
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(10a) 


[K..] = [L][L T ] 

where L is a lower triangular array. Therefore, we can write 

[K..]" 1 = [L ]" T [L]” 1 (10b) 


From Eqs . (8a) and (9), we can express the displacement in 
the modified structure at the i^ degree of freedom as 


K1 ■ l u 4 + K 


(ii) 


where the change in displacements { 6u^j is written as 


jsu*! =-[L]" T [T ]{u* 
) l J e ] e 


( 12 ) 


,-l 


and [T e ] = [L] [K ie ] . The displacement at the e^ degrees 
of freedom are obtained by solving Eq. (8b) and using substitu- 
tions from Eqs. (11) and (12), i.e. , 


[K - T T T ]{u*| = |f } - [K. ] T ju.l 
ee e e | ej l e J re j ij 


(13) 


Solving for {u e } from Eq. (13), we can determine the 
changes in the remainder of the structure by using Eq. (12). 

The above procedure has the advantage of not having to re- 
formulate the stiffness influence coefficient matrix [K ] . It 
does involve determining the value of [T e ], from Eq. (12) (a 
forward solution of the triangularized array) ; a solution for 
[u e ] (solving a system of n simultaneous linear equations, 
where n is th^ number of modified degrees of freedom) ; and 
determining (u^) from Eq . (12) (a backward solution of the 
triangularized system) . The number of computational operations 
.required for this technique, as determined in Ref. 5, is 
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( 14 ) 




where M =■ the operation count, n = the total number of degrees 
of freedom of the modified structure,. n r = degrees of freedom 
to be modified, and b = semibandwidth of stiffness matrix. The 
break-even point of this technique versus a complete reformula- 
tion occurs when n r = 0.75b. 


In Refs. 4 and 5 the condition that (F e ) =0 is imposed 
as the requirement that the "e" degrees of freedom are free to 
displace. This condition, although necessary, is not sufficient 
to determine the actual displacements associated with the modi- 
fied degrees of freedom and their influence on the remainder of 
the structure. - •• 

The algorithm for the incremental procedure to treat the 
intermittent contact and separation remains the same as that 
previously presented for the unmodified structure with the ex- 
ception that the flow of operations is interrupted between 
steps 1 and 2. Nodal displacements at nodes, chosen a priori, 
are monitored to determine whether they are to be released 
(representing crack extension), remain open (crack opening) , or 
remain fixed (crack closure) . The displacements associated with 
these decisions are then computed and their effect on the dis- 
placement vector (u} , computed in step 1, is finally deter- 
mined. 



3 . NUMERICAL RESULTS 


A uniformly loaded rectangular center-crack panel, shown in 
Fig. 1, is chosen to demonstrate the previously described method of 
substructuring to treat crack extension and intermittent contact of 
crack surfaces. The finite -element idealization of a quadrant of 
the panel, shown in Fig. 2, consists of 396 elements having 249 
nodes resulting in 474 degrees of freedom when symmetry boundary 
conditions are applied. 

The crack displacement profile, assuming elastic behavior, is 
illustrated in Fig. 3 for two crack lengths. The solid curve cor- 
responds to the displacement profile for the initial crack length 
2a; the dashed curve is associated with crack length 2a 7 , where 
a 7 /a = 1.1. These results were obtained without utilizing the sub- 
structuring technique. The circles correspond to the displacement 
profile for crack length 2a 7 , determined by means of the sub- 
structuring technique — i.e., nodes ABC and D were "broken," 
extending the crack length from 2a to 2a 7 . As can be seen from 
the figure the results from this technique are identical to those 
obtained by assuming an initial crack length of 2a 7 . 

The crack displacement profile for elastic-plastic behavior 
is shown in Fig. 4. The results are for an elastic-ideally plastic 
material and for a loading (S max ) corresponding to 31 percent of 
the yield load (Syi e id) of the uncracked specimen. At S max the 
crack is extended from a Q to a 7 (node B to C in Fig. 2) 
where a 7 /a 0 = 1.024. The results indicate a chisel-shaped profile 
for the displacements in the vicinity of the extended crack-tip. A 
similar discontinuity for the profile of an extended crack is pre- 
sented in Ref. 11. Details of the displacement profile in the 
near-tip region is presented in Fig. 5. The results for the ex- 
tended crack are compared to corresponding results for the crack 
of initial length 2a 7 . The results from the latter case do not, 
of course, indicate the discontinuous chisel-shape profile corre- 
sponding to the extended crack. 


4, CONCLUDING REMARKS 


Application of a substructuring technique to the problem of 
crack extension and closure has been outlined and implemented into 
an existing nonlinear finite element analysis program for two 
dimensional membrane stressed structures. The method, readily im- 
plemented without a significant degree of disruption of the flow 
of the original program, appears to be particularly well suited 
for adaptation within the framework of the initial strain approach 
for the treatment of nonlinear material behavior. The advantage 
is associated with the fact that the original stiffness influence 
coefficient matrix for the unmodified structure need not be al- 
tered at any point in the analysis. 

Results, demonstrating the technique to the problem of elas- 
tic and elastic-plastic crack extension, are presented for a uni- 
formly loaded center -crack panel. A more comprehensive application 
of the technique to crack extension and closure is presented in 
Ref. 12. 



APPENDIX A 


DESCRIPTION OF FAST 


FAST is an acronym for Fracture Analysis of STructures. 

The analytic foundation of the program is the displacement 
method of finite element techniques for structural analysis. 

The program represents a spin-off of a previous program labeled 
PLANE (Ref. 3) developed for the nonlinear analysis of struc- 
tures subjected to plane stress conditions. The additional 
capability of FAST, and one that distinguishes it from PLANE, 
is the ability to treat variable restraint conditions so that 
consideration can be given to the problem of an extending crack 
or to determining the effects of crack closure. 

v; 

The program is capable of treating the elastic and the 
elastic-ideally plastic response of orthotropic materials. In 
addition, consideration is given to isotropic materials exhibit 
ing elastic-ideally plastic, linear strain hardening, or nonlin 
ear strain hardening behavior. Further, the kinematic harden- 
ing theory of plasticity is used (Refs. 13 and 14) so that pro- 
vision for cyclic loading conditions involving reversed plastic 
deformation is included. 

The finite element library and a description of each ele- 
ment follows . 

Consta n t Strain Triangle 

Constant Strain Triangle (CST) , a well-known plane stress 
membrane element, was used successfully for the idealization of 
structures presented in Ref. 9. Its derivation is based upon 
the assumption of a linear distribution for the in-plane dis- 
placements u and v, and Consequently, leads to a constant 
strain state within the element (Fig. A-l) . Each vertex is al- 
lowed two degrees. of freedom (the in-plane displacements u 
arid v) for a total of six degrees of freedom for the element. 
Consistent with the total strain distribution, the initial 
strains (plastic strains) are assumed to be constant within 
each element. Stiffness and initial strain matrices have been 
developed and successfully used in Refs. 9 and 10. 



Linear Strain Triangle 

In regions of high strain gradient, the CST triangle is not 
sufficiently accurate to be used in a plasticity analysis unless 
a very fine grid is employed. The linear strain triangle (LST) 
remedies this shortcoming of the CST element. The assumption of 
a quadratic distribution for the in-plane displacements allows 
for a linear strain variation within the triangle. Two degrees 
of freedom at each node (u,v) for each of the six nodes (three 
vertex and three midside nodes) give this element a total of 12 
degrees of freedom. The initial strains are assumed to be con- 
stant within each element and are evaluated at the centroid. 

Both stiffness and initial strain matrices have been developed 
and successfully used in Ref. 10. 

Hybrid Triangles 

In transition regions — those regions in which stresses 
and strains change from rapidly varying to slowly varying — it 
becomes convenient and efficient to switch from linear strain 
triangles to constant strain triangles (Fig. A-3) , This is ac- 
complished by using four and five node triangles to maintain 
compatibility with both the CST and LST elements. These ele- 
ments together with the CST and LST elements were originally 
used in Ref. 15, and are referred to as the TRIM 3 through 
TRIM 6 family. For these mixed formulation hybrid elements, 
the displacements along edges may vary quadrat ically or lin- 
early, depending upon whether an LST or CST triangle is contig- 
uous to the respective sides. Again, the plastic strain dis- 
tribution is assumed constant within each element. 

Stringers 

For many aircraft structures — fuselages, wings, etc. — 
local stiffening is required to provide adequate stability in 
compression (Fig. A-4). Special one dimensional finite-elements 
are required to represent the stringers used for this purpose. 
Uniform cross section stringer elements have been developed 
using both constant and linearly varying strain assumptions, so 
that compatibility with both the CST and LST elements can be 
maintained. For the constant strain stringer, a linear axial 
displacement and a constant initial (plastic) strain distribu- 
tion within the element are assumed. The linear strain stringer 
stiffness matrix is based upon a quadratic displacement assump- 
tion and the initial (plastic) strains are constant within the 
element . 
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APPENDIX B 


USERS MANUAL 


The input to the program can be categorized in sections. A 
description of the sections follows. 

1 . Problem Title (20A4) 

Any 80-character title describing the problem. 

2. NPNTC NPRNT (2T5) 

0 < NPNTC <63: 

NPNTC is the sum of the following codes according to the 
options desired. A zero value of NPNTC indicates that none of 
the intermediate output will be printed. 

1 Print the load vector. 

2 Print element stiffness matrix without stringer 
contribution for elements with stringers. 

4 Print element stiffness matrix before condensa- 
tion to 8x8 or 10 x 10 for 4- and 5-node 
elements, respectively. 

8 Print each element stiffness matrix in its final 
form. 

16 Print each element stiffness matrix entry to be 
stacked with its stacking index. 

32 Print the master assembled stiffness matrix for 
the structure to be analyzed. 

NPRNT: 

If < 0, perform elastic analysis only. 

If >0, perform plastic analysis, printing output 

every NPRNT increment of load. 
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3 . Nodes Within Partitions (1615) 

A partition is defined as a set of degrees of freedom of 
a subregion of the structure that can interact only with its 
neighboring subregions. The resulting stiffness influence co- 
efficient matrix is blocked by partitions into a tridiagonal 
form as shown in Eq. (C-l) of Appendix C, where those parti- 
tions not in the tridiagonal region are null arrays. 

Specify the nodes of each partition by listing them in con- 
secutive fields. Node numbers must be positive and not exceed 
32,768. One blank field must separate the listings of each pair 
of adjacent partitions. Two consecutive blank fields (with one 
or both on a final blank card if and only if necessary) denote 
the end of the section. Two shorthands are allowed. Specifying 
m and -n consecutively is the equivalent of m, m+1, mf2, 

. . . , n, and specifying m, -n, and -p consecutively is the 
equivalent of m, mfn, nrt-2n, . . . , mfkn, where mfkn is the 
highest integer of this form less than or equal to p. 

The maximum number of partitions has, arbitrarily, been set 
to 20 and a minimum of 2 partitions is necessary for a success- 
ful solution. Consecutive numbering of nodes within any parti- 
tion is not required. It should be noted that the nodes of any 
element must be in either one partition, or at most, two con- 
tiguous partitions. The total maximum number of nodes for any 
structure is set at 900. The total number of degrees of free- 
dom (~ 2 x number of nodes) within any partition is determined 
by the value of NCORP defined in a data statement in subroutine 
MAIN. A typical value of NCORP = 10,000 words results in a 
maximum of « 80 degrees of freedom in any of two contiguous 
partitions. A more detailed discussion of the storage alloca- 
tion is presented in Appendix C. 

4. Modified Nodes and Degrees of Freedom for Elast ic 

Analysis Only: (3(215, 15 x)) MNOD, M0D0F 

The first 15 field is associated with the node to be re- 
leased. (MNOD) from a fixed boundary condition; the second, rep- 
resents the degree of freedom, u (in the global x-direction 
and/or V (in the global y-direction) . If 



MDOF =10, u is released 
= 1, V is released 

= 11, u and V are released 

TVo separate elastic solutions for displacement strains 
and stresses will be printed. The first solution is associated 
with the case where the modified nodes are assumed to be fixed. 
The second solution represents the corresponding results when 
the modified nodes are released. A blank first 15 field ends 
this section. A blank card is required when there are no modi- 
fied nodes. 

5. Nodes and Stringers for Members 2(715,2 x,3ll) 

•t 

The first 15 field gives the member number, which must be 
positive and no greater than 32,768. The next six 15 fields 
give the nodes, beginning with a major node and continuing 
around the perimeter of the member alternately major and minor. 
The absence of a minor node must be indicated by a zero or blank 
field in the proper position. The 311 digits indicate, respec- 
tively, the presence or absence of a stringer connected to sides 
4, 5, and 6, where side 4 is the side opposite the first major 
node given in the six 15 fields, side 5 is opposite the second, 
and side 6 is opposite the third (see Fig. B-l) . It is sug-' 
gested that the digits 4, 5, and 6 thus be used in the appro- 
priate positions to denote the presence of a stringer on the 
corresponding side, although any nonzero digit will suffice. A 
zero in an II field denotes the absence of a stringer from the 
corresponding side. TWo members may be specified per card 
(columns 1-40) . A zero or blank first card field ends the sec- 
tion. The program will accept a maximum number of 600 elements. 
It should be noted that the program will accept the case of 
specifying a minor node in only one of two adjacent elements 
sharing a common side. The displacement field along this side 
of the elements will be compatible only at the two major end 
nodes and not along the length of the side. 

6 . X-Coordinates of Nodes (E15.7,13I5) 

The X-coordinates of the nodes appearing in the 
are set to the value in the E15.7 field. Any number 
ation cards may be used; their first fifteen columns 
nored. A zero or blank 15 field terminates the node 


15 fields 
of continu- 
are ig- 
list for a 
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given X-coordinate. A zero or blank first 15 field (columns 
16-20) where a new X-coordinate can be specified end the sec- 
tion. Both shorthand representations of Section 3 are allowed. 
Coordinates for major nodes only are errors, but will be ig- 
nored. 

7. Y-Coordinates of the Nodes (Same as Section 6) 

8. Stringer Information (5E15 . 7 , / , (1615) ) 


To be provided only if there are stringers . The first card 
gives, in order, 

E Young 1 s modulus 

A Cross-sectional area * 

FN If SIGZS^O, FN=n, the shape parameter used 

in Ramberg-Osgood representation of stress- 
strain behavior; 

If SIGZS=0, FN=a, the slope of the linear 
strain-hardening stress-strain representa- 
tion, i.e., a = Ej/E where E T is the 
tangent modulus . 

SIGO Yield stress 

SIGZS If =/0; SIGZS = Ramberg-Osgood parameter 
a 0 . 7 

Note: If FN=0 and SIGZS=0 the material for 

the stringer element (s) is assumed to be 
elastic-ideally plastic. 

The following cards give the end points of stringers to which 
the above values are to be assigned. Every pair of nonzero in- 
tegers appearing in this input is treated as a stringer end 
point pair. Blank fields are ignored, except that a fully 
blank card ends the end point pair input. A zero E ends the 
section. If an odd number of end points are specified, problem 
execution is terminated. 
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9 . Boundary Conditions ( 14 , II , 5x , 14 15) 

The first two fields give a boundary condition specifica- 
tion, in the order u x ,Uy. 

Zero denotes a fixed degree of freedom. 

One denotes a free degree of freedom. 

Two will result in the application of a unit 
displacement . 

The 1415 fields give applicable nodes, with both shorthands of 
Section 3 permitted. Any number of continuation cards may be 
used for a given specification. A zero or blank 15 field ends 
each node listing. A zero or blank first 15 field (columns 11- 
15), where a new specification is expected, ends the section. 

If a node's boundary conditions are not specified in this sec- 
tion, both degrees of freedom are assumed to be free. 

10. Material Properties (5E15.7,/,4E15.7,/,5E15.7,/, (1615)) 

The first three cards specify material properties, as fol- 
lows . 


Card 1: 


EONE = Young's modulus in principal property axis (1) 
ETWO = Young's modulus in principal property axis (2) 
BETA = Angle measured from global x-axis and 
principal property axis (1) 

GONTO = Shear modulus in (l)-(2) principal property 
axes. 

VONTO = Poisson's ratio, v-^ 


Note (1): v^2 



Card 2: SIGOX = Yield stress in global x-axis 

SIGOY - Yield stress in global y-axis 

SIGOZ = Yield stress in global z-axis 

SIGXY = Shear yield stress in global x-y plane 


Card 3: RMOSN = If SIGZS^O; RM0SN=n, the shape parameter 

used in Ramberg -Osgood representation of 
stress -strain behavior. 
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If SIGZS=0; RMOSN=a the slope of the linear 
strain-hardening stress-strain representa- 
tion, i.e. , a = E-p/E where E>j is the 
tangent modulus. 

SIGZS = If ^0; SIGZS = Ramberg -Osgood parameters 
a 0.7 

Note (2): If RM0SN=0 and SIGZS=0 the material for the 

element (s) is assumed to be elastic-ideally 
plastic. ' 

RMOSE = Ramberg -Osgood parameter E (Young's modulus) 
YLDST = Yield stress in tension 
YLDSC = Yield stress in compression 

Note (3) : Only initially isotropic materials can be 

treated when considering linear or nonlinear 
strain hardening. 

Succeeding cards give applicable members; both shorthands of Sec 
tion 3 are permitted. Any number of continuation cards may be 
used for a given specification. A zero or blank 15 field ends 

each member listing. A zero or blank EONE ends the section. 

\ ' • 

11 . Member Thicknesses (E15. 7, / , (1615)^) 

The first card specifies the member thickness, the 1615 
cards give applicable members, as in 10. A zero or blank thick- 
ness ends the section. 

12 . Applied Load Components (3I5,4E15.7) 

Each card gives the load components applied at a member 
side, as follows: 

First 15 field: number of node (m) 

Second 15 field: number of minor node on the side, or 

v zero or blank if there is no minor node 

Third 15 field: number of other end point node (n) 

First E15.7 field: x load component at node m 

Second E15.7 field: y load component at node m 

Third E15.7 field: x load component at node n 

Fourth E15.7 field: y load component at node n 

As many load components as desired may be specified. A zero or 
blank first 15 field ends the section. 



13 . Members to Be Printed - Elastic Solution (1615) 

Specify the members whose strains and stresses are to be 
printed in the elastic analysis in the order in which they are 
desired. Blank fields are ignored, provided one nonblank field 
appears on a card. Both shorthands of Section 3 are allowed. 
Should three consecutive negative integers appear, a "1" is 
inserted before the third. A maximum of 600 members may be 
specified. Members in excess of 600 and undefined member num- 
bers are ignored. A blank card or card with only zero entries 
end the section. 

14. Nodes to Be Printed - Elastic Solution (1615) 

Up to 900 nodes whose displacements are to be printed for 
the elastic analysis; as per Section 13. o 

15. Members to Be Printed - Plastic Solution (As per Section 13) 

16 . Nodes to Be Printed - Plastic Solution (As per Section 13) 

17. Modified Nodes and Degrees of Freedom for Plastic 
Analysis Only; (3(2l5,15x)) MNOD, MDOF (As per Section 4) 

A blank card is required when there are no modified nodes. 

18. Plastic Parameters (3E15.7) 

In order, 

PPCT Load increment as a percentage of yield load. 

PMAX Maximum load to be applied. 

YEPS < 1.0; the amount by which YEPS is less than 
1.0 represents the tolerance requirements 
in determining whether a particular stress 
state satisfies the yield condition. 

19. Load Range for Problems Involving Modified Nodes 

Only: (A4, lx,El5.7) TEMP (l), PMAX 

PMAX represents the maximum or minimum value of the ap- 
plied loading for any half -cycle of load. Any number of half- 
cycle loadings may be applied. > The field for TEMP(l) must be 
blank for any specified value of PMAX and must read 


"ENDb," where b denotes a blank/ for the lost input card. 
If there are no modified nodes this section is ignored (no 
blank card is required) 


20. For Succeeding Load Cycles, One Card (I5,2E15 .7) 
Giving New NPRNT, PMAX, and PPCT 


Zero NPRNT signifies no new load cycle. To change proper- 
ties of any group of members, as given in Section 9, there may 
follow any number of cards (I5,5E15.7) giving I any member of 
the group whose properties are to be changed. 


RMOSN ” 
RMOSS 
RMOSE 
YLDST 


New Ramberg -Osgood parameters 


YLDSC _J 


1=0 indicates that the changes are complete. 


This section is ignored if the problem involves modified nodes. 

21. Each Problem's Input Must Be Ended with a Card 

Reading "ENDb, " where b Denotes a Blank, in Columns 1-4 

An additional "ENDb" is not required for those problems 
involving modified nodes . 
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APPENDIX C 


SOLUTION ALGORITHM 


The algorithm chosen for the solution of symmetric, posi- 
tive definite matrix equations in block tridiagonal form is the 
Cholesky method (referred to as Danilevskii method in Ref. 16) . 
This algorithm factors the total stiffness matrix into the 
product of a lower triangular array and its transpose and then 
solves a pair of triangular sets of equations. This factoriza- 
tion is possible only for positive definite matrices. 

Problem 


Solve AX = Y, where 

V *i B 1 

B 1 *2 B 2 

■ , A = . . : * 

b 2 ^ B 3 

B 3 A 4 


(C-l) 


is the block tridiagonal, positive definite, symmetric stiff- 
ness matrix and X and Y (representing the generalized nodal 
displacements and loads, respectively) are partitioned corre- 
spondingly as 


iT 4 

J 


fv 

CM 

X 


Y 2 


II 

. 

> 

X 3 


Y 3 

A 


J 


(C-2) 



Triangularization 

Assume A = LL^ with 


“2 L 3 

M 3 L 4 

and A partitioned as already indicated. Then 
A 1 = L 1 L 1 , 

B i = M i L i so that M i = b i l ‘i T 

P^2 ~ so that ^ 

V>2 = M 2 L 2 so that = B 2 L 2 T 

^ = M 2 M 2 + L 3 L 3 T so that L 3 L 3 - A 3 - M 2 M 2 

etc. ' 

These equations are used in turn to determine L^, Mj, L2, 
M2, etc., obtaining each diagonal block by the Cholesky algo- 
rithm and each off-diagonal block by solving triangular equa- 
tions . 
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Flow of Subroutines 



PIRATE is the user interface and supervisory routine. It 
uses the Cholesky algorithm, which factors the total stiffness 
matrix into LL^ (where L is a lower triangular matrix) and 
then solves a pair of triangular sets of equations. It is well 
known that this factorization is possible if A is positive 
definite. 


The factorization is carried out by subroutine ATTACK as 
follows: 

Subroutine FUTILE performs Cholesky factoriza- 
tions to construct the diagonal (tri- 
angular) blocks l>i, L2, . . .. 

rji 

Subroutine BOY forms the products 

, Subroutine GIRL solves equations L^M^ = b£; 

GIR.L calls TRIEQ for each column of 
B^ to get the corresponding column 
' of Mj 












PIRATE uses a large storage array; at a typical moment, 
the array holds one diagonal block of A or L and one off- 
diagonal block B or M. Diagonal blocks are placed in the 
beginning of the large array, with off-diagonal blocks posi- 
tioned at the end. This arrangement prevents storage con- 
flicts, provided that the array is large enough to hold any 
diagonal block along with either (not both) of the adjoining 
off-diagonal blocks. The symmetry of A makes it possible 
to store only about half of the diagonal blocks . 

Solution 

After factorization, the problem can be posed in the form 
LL^X = Y; PIRATE uses subroutine MORGAN to solve LZ = Y for 
Z, and then L^X = Z for X; X is then the required solu- 
tion because L^X = Z implies L(L^X) = LZ, i.e., AX = Y. 

The "forward" solution LZ = Y can be expressed as 



(C-4) 



MORGAN calls SKULL to form the products MfcZk and TRIEQ to 
solve the triangular equations for the Z^. 
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The "back" solution is so called because it obtains the 
solution elements in reverse order; in partitioned form the 
process is 


T T 

l: m. 



r Z- 

1 1 


1 


1 . 

L T m t 
L 2 2 

/ 

CM 

X 

L J 

Z 2 

T W T 

\ 


r 

/ 

L 3 m 3 


X 3 


Z 3 

T 





1 L 4 j 


^ X 4> 




(C- 6 ) 


i.e. , 


T 

L 4 X 4 


= Z, 


T 

L 3 X 3 


= Z, - 


M 3 X 4 


(C-7) 


i? 2 x 2 = Z 2 - m t 2 x 3 
etc . ' . 

X 

Once again, SKULL generates the products M^X^-n and TRIEQ 
provides the solutions X^. Since the X^ are obtained in re- 
verse order, it is necessary for MORGAN to read them backwards 
in order to produce the solution in normal order in core. 


Data Set Usage 

PIRATE uses 5 data sets, referred to here as T^, ..., T 5 . 
Ti must contain the A matrix and T 2 the Y matrix when 
PIRATE is entered. The L matrix is written on T 3 and Z on 
T 4 , with the solution X returned on T 5 in reverse order. 

It is possible to overlap the usage somewhat, but the following 
restrictions must be complied with:. 



must be different from T2 
must be different from 
must be different from 
and must be different. 


and iy 

and T, 
4 

and T,- 
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Fig. 2 IDEALIZATION OF UPPER RIGHT QUADRANT QF PANEL 





Fig 3 CRACK OPENING DISPLACEMENT PROFILE- ELASTIC BEHAVIOR 




Fig 4 CRACK OPENING DISPLACEMENT PROFILE - ELASTIC PLASTIC BEHAVIOR 
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Fig 5 COMPARISON OF CRACK OPENING DISPLACEMENT PROFILE 

ELASTIC PLASTIC BEHAVIOR 
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